The Bak-Chen-Tang Forest Fire Model Revisited 
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We reconsider a model introduced by Bak, Chen, and Tang (Phys. Rev. A 38, 364 (1988)) as 
a supposedly self-organized critical model for forest fires. We verify again that the model is not 
critical in 2 dimensions, as found also by previous authors. But we find that the model does show 
anomalous scaling (i.e., is critical in the sense of statistical mechanics) in 3 and 4 dimensions. We 
relate these results to recent claims by A. Johansen. 



During the last ten years, the concept of self-organized 
criticality (SOC), proposed by Bak, Tang and Wiesenfeld 
in Q , has been applied to a large number of phenomena. 
Although a generally accepted and rigorous definition of 
SOC still doesn't exist, a number of general features are 
supposed to hold for any model which shows SOC (see 
also §): 

(i) There should be scaling laws which in general will be 
anomalous for local models (for mean field or random 
neighbor models the exponents in general will be integer 
or half- integer, i.e. 'normal'). These scaling laws should 
not be 'trivial' (such as m ~ I 3 for the mass-length con- 
nection in 3-d), though triviality is not always easy to 
define. 

(ii) There should be no control parameter which has to 
be fine tuned. Thus the scaling should be a robust phe- 
nomenon, in contrast to standard critical phenomena. 
Again this criterion is less clear-cut then one might wish. 
In particular, 

(iii) SOC shows up usually in slowly driven systems, when 
the driving rate tends to zero (an advocatus diaboli who 
insists that this rate should be considered as a control 
parameter could thus conclude that SOC doesn't exist 
at all). Typically, such systems become locally unstable 
when the stress exerted by the driving exceeds some limit, 
and react with 'avalanches' of activity which are large on 
microscopic scales, but small on macroscopic ones. 

(iv) Finally, a common feature of many systems with 
SOC is that the driving is not controlled as usual through 
a force, but through a flux (i.e., an extensive quantity) 

Phenomena where these features show up include sand 
piles , earth quakes H|| , pinned surfaces J|,[?J , and bi- 
ological evolution ||. A last application are forest fires 
with small growth rate of trees and even smaller rate for 
spontaneous ignition ('lightning') ||. The latter example 
is special in the sense that it requires three different time 
scales for criticality. 

In |M (BCT), it was claimed that this list should also 
include forest fires without lightning. The specific model 
studied by BCT used a regular <i-dimensional lattice and 
discrete time, with each lattice site in one of 3 possible 
states: green tree, burning tree, or ash. During one time 



unit, a burning tree ignites all green neighbors (if there 
are any) and turns itself into ash. This is the fast part 
of the dynamics. The slow part is the re-growth of trees. 
It is modelled by a stochastic spontaneous transition ash 
— > tree with probability p< 1. Thus in each time step a 
randomly selected fraction p of all ash sites is flipped into 
trees. As pointed out in [jn), a much better interpreta- 
tion of the BCT model would be in terms of epidemics 
with slow recovery resp. slow loss of immunization, but 
we shall continue to speak about forest fires for conve- 
nience. 

Simulations jyj on large lattices (up to 4800 2 sites) 
and for very small values of p (down to 5 x 10~ 4 ) showed 
however that dynamics in d = 2 is quite different from 
that suggested by BCT on the basis of small-scale sim- 
ulations. Instead of being critical, the system develops 
noisy spiral patterns which become less and less noisy 
as p — ► 0. More precisely, spiral arms (fire fronts) propa- 
gate with finite mean velocity v for any p, and the typical 
distance £ between spiral arms (the characteristic length 
scale) as well as the time T between two passings of a 
front scale as 1/p: £ cx T oc 1/p, In the limit p — ► the 
coherence length thus diverges, implying that the dynam- 
ics is governed by average tree densities over larger and 
larger regions. In this limit the dynamics is thus similar 
to that of coupled relaxation oscillators with very sudden 
discharge and very slow recharging. This explains qual- 
itatively the patterns found, though any details are still 
badly understood due to the inherent noise for any p 
which leads to permanently ongoing pattern rearrange- 
ments. 

Indeed, the authors of |ll| were rather careful in the in- 
terpretation of their simulations. They pointed out that 
if the model were non-trivially critical with characteristic 
length and time scales scaling as 



£ oc p a , T cx p 



-0 



a, 0^1, 



(1) 



then this could only be the case if fire fronts get fuzzy for 
p — ► 0, and front velocities would tend to zero. In such a 
case the fire would be at a critical percolation threshold 
(a similar scenario does hold indeed for some versions of 
forest fire models with lightning [ p^ , p"3| ) . The medium 
into which it percolates would not be uncorrelated, in 
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contrast to standard percolation. Using nevertheless es- 
timates of critical exponents from the latter, it was con- 
cluded in jflj that such an alternative scenario was very 
unlikely but hard to exclude rigorously. The same con- 
clusion (with less caveats in spite of less statistics) was 
reached also in [H . 

Simulations in 3 and 4 dimensions were much less sig- 
nificant. This was partly due to difficulties in visualizing 
such systems, partly because it is hard to keep a fire from 
getting extinct for small values of p. Nevertheless some 
indications for well defined fire fronts were seen in d = 3 
p| , and it was concluded that basically the same sce- 
nario holds there as in d = 2. Obviously, this cannot 
be true for arbitrarily high dimension. At least for d > 6 
the behavior should be that of dynamic percolation, with 
small fluctuations of the density of trees around a mean 
field value. The characteristic time of these oscillations 
is not T ~ l/p as for fronts, but T ~ l/^/P 101 ■ Simula- 
tions showed the actual behavior for d = 3 and d = 4 was 
in between both ]ll| , |T4| . Unfortunately, this was not fol- 
lowed up, and the possibility T ~ l/p 13 with 1/2 < < 1 
was not considered seriously. 

In a recent paper, Johansen claimed exactly that. 
In addition, he claimed that the same is true also in d = 2. 
On the other hand he confirmed that spirals are formed 
in d = 2, and that the typical distance between fire fronts 
scales as L ~ l/p. Now it is easy to see that the latter 
statements are self-contradictory. They would imply that 
fronts propagate with speed v = L/T ~ p-i 1 -/ 3 ) — > oo 
for p — > 0. Since the front can propagate at most one 
lattice site in each time step, this is impossible. 

In order to clarify this situation, we report in the 
present letter on simulations where we measured T with 
high precision, for d — 2,3 and 4, and for wide ranges 
of p. For d = 3 we find indeed a very clear indication 
of anomalous scaling, (3 — 0.77 ± 0.02. The situation is 
slightly less clear in d = 4 for reasons detailed below, 
but we again find scaling (with (3 « 0.6 ± 0.05). These 
findings imply that pT — ► for p — > 0. For d = 2, finally, 
we find that pT also decreases slightly with decreasing p, 
but not as a power law. Our data are not precise enough 
to distinguish clearly between a logarithmic increase, 

pT~[l/log(l/p)P, 7 >0, (2) 

and a limited increase which leads to a finite value at 
p — > 0. They favor the latter. But if eq.(||) would hold, 
instead, we would have the alternative scenario men- 
tioned above in which fronts become fuzzy, front veloci- 
ties become zero, and the evolution is basically a critical 
(correlated) percolation phenomenon. 

The latter seems to apply in > 3 dimensions, although 
the situation is not clear there either. The problem is 
the following: if pT — > for p — > 0, then the fraction of 
replenished trees between two peaks of the local fire in- 
tensity has to go to zero also. Except for transients this 



means that also the fraction of trees burnt during such a 
peak must tend to zero. In such a case we would naively 
expect that the amplitude of any (noisily) periodic ob- 
servable should diminish when p is decreased. This is not 
observed. Instead, the peaks stand out very clearly, even 
for the smallest values of p. Thus our data suggest at 
first sight that fires burn large areas completely, also in 
3 and 4 dimensions (they do so in d = 2, but there it is 
expected). This might indicate that even with p values 
as small as 10~ 3 we are not yet in the asymptotic scal- 
ing region, but an alternative scenario will be discussed 
below. 

For our simulations we used the basic algorithm de- 
scribed in jfl). The current state of the system is en- 
coded in two data structures: a list of burning tree sites, 
and a bit map indicating for each site whether it contains 
a tree or not. Notice that we do not have to distinguish in 
the latter between burning and non-burning trees since 
that information is contained in the list. Thus we can 
use bit coding in order to simulate very large lattices. 
Instead of using d coordinates, each site is indexed by a 
single integer, and boundary conditions are helical. In 
each time step, a new list of burning trees is established 
by scanning through all neighbors of all entries in the old 
list, and the old list is thereafter replaced by the new one. 
After this, pN^ sites (iV as h is the total number of sites 
containing ash) are randomly selected and switched from 
ash to tree. 

To avoid that the fire dies out, we used correlated ran- 
dom initial conditions, and discarded transients which 
involved at least 100 oscillation periods. If the fire died 
out nevertheless, we ignited some fires 'by hand', and 
started the run again. 

Lattice sizes went up to 16384 x 16384 in 2 dimensions 
(for p = 0.00005), and to comparable numbers of sites in 
d > 3. Simulation times went up to t — 4 x 10 6 , and were 
always larger that 100/p. 

Since the aim of the present paper was to obtain pre- 
cise estimates of T, the biggest effort was devoted to that. 
As in previous analyses, we used the number of burning 
trees as observable. But in order to improve the signal 
to noise ratio, we did not simply measure the total num- 
ber, averaged over the entire lattice. The reason is that 
fires in different parts of a large lattice will in general not 
oscillate in phase, whence cross correlations from distant 
regions will mainly contribute to the noise in the auto- 
correlation function. We therefore proceeded as follows: 
we divided the lattice into hypercubes of linear size I with 
/ < L, and measured the numbers rij(t) of burning trees 
in the ith cube at time t. From each of these local time 
series we estimated autocovariances 

Ci(jt) = (m{T)ni{t + T)) (3) 

which were then averaged over the lattice, 
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FIG. 2: Log-log plot of pT versus p for 2 dimensions. Here, 
T is the average peak-to-peak distance in c(t). 



Oscillation periods T were either estimated from peak- 
to-peak distances in c(t) or by Fourier transforming c(t) , 
obtaining thereby maximum entropy spectrum estimate 
S(f) |l6|| . Results were the same. Two typical plots of 
c(t) are shown in fig. 1. 
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FIG. 1: Autocovariances c(t) with: (a) d = 2, p = 0.001, 
lattice size 4096 x 4096, 10 6 iterations, and I = 128; (b) d = 3, 
p = 0.003, lattice size 256 3 , 10 s iterations, and I = 16. In both 
cases, only 16 squares resp. cubes were used in the averaging. 
Normalization is arbitrary. 

The sharpness of the peaks in c(t) and the strong 
higher harmonics result from the fact that I << £, 
whence rij(i) is non-zero only during the short time when 
a fire front passes through cube i. As a consequence we 
obtain very clean signals and very precise estimates of T. 
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We should point out here that T is not the only time 
scale characterizing c(t). First of all, there is also the co- 
herence time T co h- It can be measured from the asymp- 
totic exponential decay of the oscillation amplitudes. As 
expected, it grows quickly with I /p. But we did not make 
systematic measurements since the exponential decay is 
not observed at finite times due to a third time scale, 
namely the regeneration time r reg cn = It is easy to 
see that the amplitudes of all peaks at t > have to de- 
crease for fixed I and p — *■ 0, if T <C T regen in this limit. 
Indeed, for I = 1 one finds c(T)/c(0) « pT. Thus the 
shape of c(t) depends strongly on the box size I, but we 
verified that the locations of the maxima (which deter- 
mine T) are independent of it. 

For small values of p (requiring large systems and long 
simulation times) it was not feasible to store all rii(t) for 
each cube i and every t. In such cases we decimated the 
data by either reading out only a fraction of cubes, or by 
coarse-graining in t and storing n'^kt) — Yl^kt* n i( T ) 
for some integer k > 1. Both gave nearly the same per- 
formance as without decimation, even if the data were 
reduced by more than one order of magnitude. 
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FIG. 3: Same as fig. 2, but for d = 3. The dashed line 
corresponds to a power law T ~ p~^ with f3 — 0.77. 

Our final results are shown in figs. 2 to 4. Each of them 
is a log-log plot showing pT versus p. In figures 2 and 3 
the estimated errors are smaller than the symbols. We 
see clearly the trends discussed above. For d = 4 (fig. 4) 
the errors are larger and not purely statistical: runs with 
different initial conditions gave occasionally values which 
differed by more than naive statistical error estimates. 
Obviously this means that either the system is not er- 
godic, or that our simulation times were not sufficient to 
explore all phase space. Nevertheless we believe that the 
data give a clear indication of scaling also for d = 4. 
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FIG. 4: Same as fig. 2, but for d = 4. The dashed line 
corresponds to /3 = 0.6. 

As pointed out above, the most straightforward inter- 
pretation of our data for d > 3 is that there are fronts 
with distances £ ss vT ~ Between two successive 

fronts only a tiny fraction of trees can re-grow, whence 
also a tiny fraction of sites could burn when the front 
passes. If these fronts would propagate into an essen- 
tially unstructured medium, this would imply that the 
process is close to critical percolation, and it would be 
hard to understand the very large and regular amplitudes 
seen e.g. in fig. lb. To understand better what is going 
on we tried to visualize typical 3-d configurations, but 
only with moderate success. But the above suggests that 
fronts do not propagate into an unstructured medium. It 
is well known that the complement of a slightly super- 
critical percolation cluster in d > 3 is connected. Thus 
a passing fire, even if it is supercritical and has thus a 
sharp front, could still leave intact connected regions with 
slightly subcritical or even supercritical tree densitites. 
This would allow the next front to pass very soon again. 
The crucial point is that this front should propagate only 
on a sparse but connected subset of trees, requiring the 
patters of trees to be highly structured. 

In summary, we have shown that there are anomalous 
scaling laws in the Bak-Chen-Tang forest fire models, but 
only in higher dimensions. This suggests that there are 
sharp fire fronts in d > 2, even in the limit p — ► 0, but 
each front propagates only on a sparse subset of trees. In 
this way the fire can be endemic in d > 3, burning only 
an infinitesimal fraction of trees between two recurrences 
to the same mesoscopic region without, nevertheless, re- 
sembling critical percolation. 

Since the present model can be considered as a model 
for extremely noisy coupled relaxation oscillators, it is an 
interesting question whether this is the generic behavior 
of noisy coupled relaxation oscillators in > 2 dimensions. 

Another final remark concerns the relationship of the 



present model with other SOC models. The main dif- 
ference is that in the present model the system is slowly 
driven (by the growth of trees) into a susceptible state, 
while most other SOC models are driven into unstable 
states which will 'topple' (discharge, catch fire, ...) spon- 
taneously. This difference would be blurred if a suscep- 
tible site has a small chance to topple anyhow, as in the 
Drossel-Schwabl H model. But it is also blurred if the 
connectivity of the lattice is so high that activity can effi- 
ciently spread over large distances without leaving many 
traces. This is obviously what happens in the present 
model when d > 3. This explains why the model shows 
SOC for d > 3, but not for d = 2. 
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